Molecular simulation analysis of structural variations in lipoplexes 
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Abstract 

We use a coarse-grained molecular model to study the self-assembly process of com- 
plexes of cationic and neutral lipids with DNA molecules ("lipoplexes") - a promising 
nonviral carrier of DNA for gene therapy. We identify the resulting structures through 
direct visualization of the molecular arrangements and through calculations of the cor- 
responding scattering plots. The latter approach provides a means for comparison with 
published data from X-ray scattering experiments. Consistent with experimental results, 
we find that upon increasing the stiffness of the lipid material, the system tends to form 
lamellar structures. Two characteristic distances can be extracted from the scattering plots 
of lamellar complexes - the lamellar (interlayer) spacing and the DNA-spacing within each 
layer. We find a remarkable agreement between the computed values of these two quan- 
tities and the experimental data [J. O. Radler, I. Koltover, T. Salditt and C. R. Safinya, 
Science, 1997, 275 810-814] over the entire range of mole fractions of charged lipids (CLs) 
studied experimentally. A visual inspection of the simulated systems reveals that, for very 
high fractions of CLs, disordered structures consisting of DNA molecules bound to small 
membrane fragments are spontaneously formed. The diffraction plots of these non-lamellar 
disordered complexes appear very similar to that of the lamellar structure, which makes the 
interpretation of the X-ray data ambiguous. The loss of lamellar order may be the origin of 
the observed increase in the efficiency of lipoplexes as gene delivery vectors at high charge 
densities. 
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When DNA molecules are mixed with neutral and cationic lipids (CLs) in an aqueous envi- 
ronment, they spontaneously aggregate to form macromolecular complexes called "lipoplexes" . 
These complexes have attracted much attention over the past two decades because of their po- 
tential use as nonviral transfection vectors in gene therapy [TH5]. Transfection is a two-stage 
process involving adsorption and entry (via endocytosis) of the lipoplex into the cell, followed 
by the release of the DNA to the cytoplasm and delivery to the nucleus, which makes the DNA 
available for expression [BHH]. CL-DNA complexes exhibit low toxicity and nonimmunogenicity, 
but their transfection efficiency (TE) remains low compared to that of viral vectors [HE]- This 
has spurred an intense research activity aimed at enhancing TE. Recognizing that the struc- 
ture of CL-DNA complexes may strongly influence their function and TE, much of the effort in 
theoretical and experimental studies has been devoted to understanding the mechanisms gov- 
erning complex formation, structure, and phase behavior. X-ray diffraction experiments have 
revealed that CL-DNA complexes exist in a variety of mesoscopic structures. These structures 
include: (i) a multilamellar phase where DNA monolayers are intercalated between lipid bilayers 
(-^a) [12]' an d (ii) the inverted hexagonal phase with DNA encapsulated within monolayers tubes 
and arranged on a two-dimensional hexagonal lattice {Hfj) 

One of the major advantages of lipoplexes over viral capsids is their ease of preparation 
and their almost unlimited DNA-carrying capacity, which stem from the fact that the vector is 
formed by spontaneous self-assembly in aqueous solutions. The electrostatic attraction between 
the anionic DNA and the CLs along with the entropic gain associated with the release of tightly 
bound counterions from the CLs and DNA are the driving forces for the formation of a complex. 
In a recent publication, we reported on coarse-grained (CG) simulations of self-assembly of CL- 
DNA complexes [12] . We demonstrated, in agreement with previous theoretical studies and X-ray 
scattering experiments [10 [ll llfT3|. ll4j . that rigid membranes tend to form lamellar complexes. For 
soft membranes, the preferred geometry is that of the inverted hexagonal phase. Our simulations 
also revealed that the phase diagram of the CL-DNA complexes is quite rich and includes, in 
addition to the lamellar and inverted hexagonal complexes, several other disordered structures 
with distinct configurational characteristics. We also found a new ordered phase, which has 
thus far not been observed experimentally, where DNA rods and cylindrical micelles form a 2D 
square lattice analogous to the 3D cubic NaCl-type structure. Our analysis of the computed 
self-assembled structures was based on simulation images and on the calculation of the Fourier 
transform of the DNA positions [12]. The Fourier transform provides a quantitative measure for 
how the simulated structures would appear in x-ray scattering experiments. To achieve a better 
comparison with the experimental data, one should also consider the contribution of the lipids 
to the scattered intensity. By plotting the separate contributions of each component (something 
which cannot be done experimentally), one can dissect the information displayed in the scattering 
plots. 

In this paper we analyze the scattering intensity plots (the square of the Fourier transform 
averaged of all angles [12]) of the lamellar complexes, which are formed in our simulations 
when the DNAs are mixed with stiff lipid material. The technical details of the model and the 
simulations have been presented in ref. [12]. In short, the model is based on the Noguchi-Takasu 
implicit solvent CG membrane model [T5] in which each lipid of length /lip is represented by 
a linear rigid molecule [16J consisting of three beads of diameter a = /lip/3 = 6.25A, one of 
which is hydrophilic and the other two are hydrophobic. The CLs are modeled by associating 
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the hydrophilic bead with a positive unit point charge, while the DNA molecules are modeled 
as infinitely long parallel rigid rods of diameter -Ddna = 4<t = 25A with uniform charge density 
corresponding to -1.7e/A. The molecules interact via three types of interactions: (i) Unscreened 
electrostatic interactions which are calculated using the Lekner summation method [TTIHH] - (ii) 
Short-range repulsive ("hard core") potential. The bead-bead pair potential U rep ^ is given by 
Eq. (4) in ref. [To], and the bead-DNA potential £7 r ep,bD( r ) = ^re P ,bb( r — 1.5er). Since the DNA 
rods are strongly repelled from each other by electrostatic forces, there was no need to introduce 
an additional £/ r ep,DD- (hi) The Noguchi-Takasu hydrophobic interaction potential, given by 
Eqs. (5)-(6) in ref. [15] . The lipids and the DNA rods are initially randomly distributed within 
a given volume in the simulation box and, through Molecular Dynamics (MD) simulations at 
constant temperature, we follow the evolution of the complexes under different conditions. As the 
objective of the simulations is to attain self-assembled structures representative of equilibrium, 
we simulate tens of millions of time steps such that at least half of the total simulated time does 
not change the characteristics of the visible structures. We have also verified that, while the 
details of the shown structures do depend on initial conditions, the significant characteristics, 
such as the peaks in the resulting scattering intensities, are well defined. 

We study isoelectric complexes where the total charges on DNA and the CLs neutralize each 
other, with no added counterions. The structure of the complex is determined as a function of 
two parameters: (i) the fraction of CLs, C , which can be varied by adding different amounts 
of neutral lipids (NLs), and (ii) the bending modulus, n s , which is the prefactor of the bending 
energy term introduced in the later version of the Noguchi-Takasu model [T9] to control the 
stiffness of the simulated membranes (see Eq. (7) in ref. [12]). 

Fig. [1] shows the diffraction patterns of stiff complexes (/t s = 10) with different values of 
<f) c ranging from C = 1 (a) to C = 4/15 (1). Following the approach outlined in ref. [12], we 
calculate diffraction patterns from the two-dimensional (2D) Fourier transformation 

N 

= ^2 w j exp(ifj-g) , (1) 

where fj represents the 2D coordinate of the DNA rod or the 2D coordinate of a bead's center of 
mass in the plane perpendicular to the DNA axis, q is the reciprocal vector, and Wj represents 
the electron density of the jth particle relative to bulk water. For each value of <p c , we show a 
triplet of figures consisting of (left) the simulated scattered intensity from the DNA rods (with 
iV = Ndna in Eq. ([I])), (middle) the simulated scattered intensity from the lipids (N = 3A^lip), 
and (right) the self-assembled structure of the complex. The displayed scattering intensities 
are I(q) ~ ((^-'(g)! 2 )^, where q = qexp(i9) and (. . ) e denote the average over all angels. All 
the DNA scattering plots are drawn on the same scale. The lipid scattering plots are also 
drawn on the same scale, except for (a)-(e) which are multiplied by the factor indicated on the 
corresponding plot. The relative scattering intensity of the lipids and the DNA depends on their 
electron densities, as well as on C . Using reasonable values for the electron densities [20|[2Tj . 
we find that the scale of the DNA intensities is two order of magnitude larger than that of the 
lipids. Nevertheless, the scattering plots of the lipids exhibit two well identified peaks located 
at q^AM and 2 (/lam- These peaks are commonly associated with the lamellar structure. From 
the position of the first lamellar peak, one can extract the inter-layer lamellar spacing d through 
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Qlam = 2n/d. The DNA plots generally exhibit three peaks, two of which coincide with the 
lamellar peaks from the lipid scattering plots, and one which is commonly referred to as the 
"DNA peak". The position of the latter at <7dna provides information about the DNA spacing 
within each layer of the lamellar structure, g?dna- ft is assumed that ^dna = 2%/ g?dna- Notice 
that the scattering from the DNAs includes the information about the lamellar spacing observed 
also in the scattering intensity from the lipids. This can be easily understood by considering the 
idealized lamellar structure sketched in fig. [2j In this structure, the DNA rods form an oblique 
lattice with lattice vectors equal to a\ = g?dna and a-i = dj sin#. The reciprocal lattice is also 
oblique with the same angle 9 between the lattice vectors b% = 2ir/d and 62 = 27r / g?dna- The 
peaks which can be seen in the scattering plots of the DNAs correspond to the reciprocal lattice 
vectors: bi, 2bi, and 62- These are also the peaks which are usually observed in actual scattering 
experiments [10]. For idealized lamellar structures, the scattering intensity has peaks at other 
wavevectors q which correspond to linear combinations of integer- multiples of &i and 62- The 
positions of these peaks are ^-dependent [22]. In non-idealized complexes, like the ones in our 
simulations, these peaks are usually very small and are hard to be detected in the scattering 
plots. 

Two open arrows are drawn in each of the lipid scattering plots included in fig. HJ The first 
one indicates the position of the larger lamellar peak which is located at (/lam = 2n/d. The 
second one denotes the wavevector 2<2lam, where we indeed find the second lamellar peak. These 
two arrows are copied into the corresponding DNA scattering plot, verifying that these peaks 
are also reproduced by the DNA ordering as explained above. Notice that (/lam is only weakly 
dependent on <f> c . In contrast, the position of the third peak, which is indicated by the solid 
arrow, varies noticeably with C . This peak is located at (?dna — 27r/<iL>NA> and the variations 
in its position reflect the decrease in the DNA spacing with increasing C . From the computed 
scattering plots shown in fig. [TJ we can extract d and g?dna as a function of <p c . Our results are 
summarized in fig. [3] (a). The visual impression of this plot is that the results are in overall 
agreement with the idealized geometry of a lamellar structure shown in fig. [2j Specifically, the 
lamellar spacing is well approximated by the sum of twice the length of the lipids and the diameter 
of the DNA rods, d ~ (2Zlip + -Ddna) ~ 10cr = 62. 5A. Similarly, the characteristic distance c^dna 
between DNA rods in each layer approximates a linear relationship with l/<j) c [T( H l2"3 | l2"4"]. This 
relationship can be derived from the simple geometric consideration that equally spaced DNA 
rods fill the surface area made available by the lipid bilayer material for any given charge density. 
There are, however, noticeable deviations from this idealized picture, especially at large <p c . The 
lamellar spacing decreases below the ideal value, while the DNA spacing attains values higher 
than predicted by the linear relationship. The origin of these discrepancies becomes clear by 
inspection of the corresponding structures shown in fig. [TJ where we observe that the long range 
lamellar order is lost at high C in favor of a disordered arrangement of smaller DNA-bilayer 
fragments. The transition is consistent with previously described membrane rupture at high 
charge densities resulting from the electrostatic stresses that develop in the complex [23||24"]. 
The highly charged cationic membrane fragments strongly associate with the DNA rods to form 
the disordered structures seen in fig. [TJ (a)-(d). We should not, therefore, be surprised that 
the values of d and g^dna, which we inferred for lamellar structures, deviate from the expected 
behavior. The nearly constant behavior of g?dna vs. c in the disordered phase can be well 
observed in the structures seen in figjl] (a)-(d). The decrease in d in this regime is also consistent 
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with the transition into the disordered phase, where the highly charged cationic membranes 
become squeezed between the negatively charged DNA rods. 

Fig. [3] (b) shows the synchrotron X-ray scattering data reported in ref . [10] . The agreement 
with the simulation results in (a) is obvious, lending credibility to our CG model as well as to the 
Fourier space analysis of the resulting structures. The horizontal axes of the figures express the 
inverse of the membrane charge density using different scales: 1/0 C in (a) and L/D (the mass ratio 
between the lipid and DNA material) in (b). The scales are linearly related by: 2.2(1/0 C ) = L/D. 
In both figures, the lamellar spacing approaches the asymptotic value d ~ 2/ L ip + -Ddna at small 
membrane change densities. These asymptotic values are different in the two figures, but this is 
merely a consequence of the chosen model parameters for /lip and -Ddna, which slightly differ 
from the experimental values. Both figures exhibit a weak, and very similar, monotonic decrease 
of d with increasing membrane charge density. For fully charge membranes, the value of d is 
depressed by about 15-20% compared to the low density asymptote. In ref. [25] this decrease has 
been attributed to the difference in length between DOPC (neutral) and DOTAP (cationic), the 
latter being about 6A shorter than the former. We, however, observe the same behavior with 
CLs and NLs being geometrically identical. Our simulations point to two more explanations 
for this observation: For moderately charged membranes the effective bilayer thickness slightly 
shrinks with C due to the tensile stress induced by the (negative) electrostatic energy density 
of a confined charge neutral systems [17] . For high membrane charge density, the decrease in d 
is likely related to the loss of lamellar order. In this regime, the derived value of d = 2-7r/gLAM 
does not necessarily coincide with the actual interlayer spacing, since the derivation is based on 
the presumption that the complex is ideally lamellar. 

The agreement between the simulation and experimental results for g?e>na i s also clear. In 
both fig. [3] (a) and (b), we observe that the DNA spacing drops from g?dna ~ 60A at low charge 
densities to g?e>na ~ 30A at high charge densities, in a manner which is well approximated by a 
linear relationship with the inverse charge density. At very high charge densities, both figures 
exhibit the same deviation from a linear relationship between g?dna and the inverse charge density. 
This feature has been attributed in ref. [10J to the limiting contact distance, -Ddna, between DNA 
rods. This interpretation of the results is correct provided that the hydration shell is included in 
the contact distance. Our results provide yet another possibility. Visual inspection of the self- 
assembled structures show that the DNA rods do not experience any hard core interactions. The 
plateau-like behavior of g?dna at high charge densities is related to the formation of disordered 
structures which enable a more loose packing of the DNA rods. 

What we have described in this paper is based on the very close agreement between the 
computational and experimental results shown in fig. [31 These results can be explained by 
a structural shift from lamellar to fragmented geometry occurring at high membrane charge 
densities. The fragmentation of the membranes is consistent with the membrane rupture observed 
in our earlier work described in refs. [23l[2l], where we used a different CG membrane model. 
This consistency gives us confidence that the loss of structural integrity of the membrane at high 
charge densities is not an artifact of a particular model. However, we recognize that although this 
work features the largest complexes ever simulated, there may still be some finite size effects that 
obscure the comparison with experiments. One such finite size effect is related to the periodic 
boundary conditions along the DNA axis, which enforce the infinite DNA rods to lie parallel 
to each other. This constraint, which simplifies the computational scheme, may lead to the 
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formation of structures with artificial spatial correlations between the DNA rods. Another finite 
size effect is related to the relatively small sizes of the simulated complexes which, therefore, have 
scattering plots with peaks that are broader than in the corresponding experimental scattering 
plots. This low resolution makes it difficult to infer the degree of order from the width of the 
peaks. Yet another consequence of finite sizes is the enhancement of surface effects. While it 
seems plausible that the increase in the electrostatic tensile stress at high charge densities does 
proliferate structural defects, it can be that these defects form more easily on the boundaries of 
the complex and, therefore, they become over-expressed in our smaller complexes. Nevertheless, 
the close agreement between fig. [3] (a) and (b), over the entire range of charge densities, supports 
the possibility that structural variations observed in our simulations may take place in nature. 
Such a structural shift from lamellar to fragmented geometry should have implications for gene 
therapy. The shift may explain the improvement in transfection efficiency (TE) exhibited by 
these complexes at high membrane charge densities [6]. One of the main limiting stages in the 
transfection process is the release of the genetic material from the complex into the cytoplasm of 
the host cell. It is indeed reasonable to expect that the DNA rods will be more readily released 
from the fragmented disordered complexes than from lamellar structures with long range order. 
Given the remarkable success of our model, which is based on a highly CG representation of the 
constituting molecular species and their interactions, it is fair to anticipate the application of 
this model for obtaining direct observations of the mechanisms governing transfection and gene 
delivery. 
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Figure 1: Self- assembled complex structures consisting of 32 DNA rods mixed with 800 CLs. 
The amount of NLs varies from (0 C = 1) at (a) to 2200 (0 C = 4/15) at (1). Each structure is 
initiated in a random molecular configuration, and has evolved for 10 — 50 x 10 6 MD time steps. 
The structures are viewed along the DNA axes. Color coding: grey - hydrophobic lipid beads, 
red - charged hydrophilic heads, green - neutral hydrophilic heads, and blue - DNA rods. For 
each configurations, the scattering intensities of the DNA rods and the lipids are also plotted. 
The open arrows indicate the position of the lipid peaks at ^lam and 2g L AM- The solid arrow 
indicate the position of the DNA in-plane correlation peak at ^dna- 
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Figure 2: Sketch of an idealized lamellar complex. 
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Figure 3: (a) The inter-layer lamellar spacing d and the DNA-spacing g?dna vs. l/<f> c , as computed 
from the peaks indicated by arrows in the scattering plots shown in fig. [TJ (b) The same quantities 
derived from synchrotron X-ray scattering data reported in ref. [10]. 1/0 C = 1 in (a) corresponds 
to L/D = 2.2 in (b). (b) has been adopted from [T0j[25] with permission. 
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